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Functional data analysis in a mixed-effects model framework is done using operator calculus. 
In this approach the functional parameters are treated as serially correlated effects giving an 
alternative to the penalized likelihood approach, where the functional parameters are treated as 
fixed effects. Operator approximations for the necessary matrix computations are proposed, and 
semi-explicit and numerically stable formulae of linear computational complexity are derived for 
likelihood analysis. The operator approach renders the usage of a functional basis unnecessary 
and clarifies the role of the boundary conditions. 
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1. Introduction 

The aim of this paper is to derive operator approximations of the matrix computations 
used to estimate the fixed and the random effects in a mixed-effects model, where M 
samples j/i, . . . , j/a/ G of temporal curves have been observed at N predefined time 
points ti,...,tN. The main technical contribution of this paper, making it practically 
possible to solve the estimation problem as a functional estimation problem, is that the 
proposed operator approximations have linear computational complexity in the sample 
length N . Consequently, the mixed-effects inference becomes feasible in the realm of 
functional data analysis, where N can be large. 

Concatenating the samples ym = {Vmn}n=i,...,N G into an observation vector y = 
{ym]m=i,...,M G R^'"'"' with dimension iVtotai ~ N * M the statistical model we use is 
given by 

y^ri3 + Zu + x + e. (1) 

In this linear mixed-effects model the design matrices F e K^totaixp g^j^^ ^ ^ jjJVtotaix? ^re 
known and assumed to have full ranks p and q, respectively, and the fixed effects /3 € W 
and the random effects u ^ Afq{0,a^G) may be shared by the M samples. The random 
component x = {xm}m=i,...,M ^ -^Ntot^iiO, o''^ R) is partitioned in the same way as the 
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observation vector y and consists of discretized readings Xm = {a;^*(tri)}ri=i,...,Af S of 
unobserved (latent) random functions x^* : [a, b] — >■ R. We assume that the random func- 
tions x^i^ , . . . ,x^^j arc independent and identically distributed Gaussian processes with 
zero mean. The covariance matrix a^R will be specified below appealing to the smoothing 
splines methodology often used in functional data analysis. Due to the i.i.d. assumption 
there exists a covariance matrix Rq S W^^^ such that R = Rq® Im, where (8) is the Kro- 
necker tensor product, and 1^/ S ^mxm identity matrix of dimension M. The last 

component in the mixed-effects model is the measurement noise e ''^ J^Ntots^ii^T '^'^^Ntot^i) ■ 

Our objective is to derive computationally efficient formulae for the maximum likeli- 
hood estimate of the fixed effects /?, the best linear unbiased predictions of the random 
effects u gM."^ and of the latent random functions x^' : [a, 6] M, and for the restricted 
likelihood function. The latter allow for restricted likelihood inference on the variance 
parameters <t^ > 0, G gM.'^^'^ and Rq G R^^^ . The methodology presented in this paper 
has two notable differences as compared to the penalized likelihood approach to func- 
tional data analysis; see, for example, the books by Ramsay and Silverman [11, 12]. From 
the viewpoint of computations we devise methods that work directly on the data vector 
y and; for example, provide predictions i9^E[a:^';'(t)|y] of the temporal derivatives of the 
latent functional parameters. In particular, there is no basis representation of the func- 
tional object E[a:^*|y]. This is by contrast with the standard technology used in functional 
data analysis, where functional parameters are given a finite dimensional representation, 
for example, in a spline basis, and the sparseness of the associated covariance matrices 
is invoked to achieve feasible computations. As an alternative to this we use analytically 
tractable operator approximations of the matrix equations. From the viewpoint of statis- 
tical modeling we model the functional parameters x^^ as random effects. Whether this is 
preferable over the fixed effect interpretation underlying the penalized likelihood depends 
on the particular application at hand. The distinction between random and fixed effects 
is here the same as for classical mixed-effects models; see [13] for a thorough discussion 
of the issue and [6] for a comparison of the associated inference methodologies. 

In the simplified version y = x + e of model equation (1), the sample ym may be un- 
derstood as a noisy observation of the function a;^* : [a, 6] — >■ R taken at the sample points 
ti, . . . ,ti<! ■ In the penalized likelihood approach to functional data analysis the functional 
parameters x^^ are treated as fixed effects. The penalized negative log likelihood is given 



where is a differential operator of some order k measuring the roughness of a function 
9 £ C'^{[a, 6]; R). The so-called smoothing parameter A > quantifies the trade-off between 
a close fit of the observations and the roughness of the functional parameters. Since the 
space of functions is infinite-dimensional, such a trade-off is required to avoid overfitting 
of the finite number of data points. 

In this paper we avoid the curse of dimensionality by providing the theoretical solution 
in the function space before plugging in the observed grid readings to compute the so- 
lution. This is done using the operator ^ = J(^^J(^, which is of order 2k and defined on 
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C'^''{[a,b];R). To ensure positive definiteness of ^ we impose boundary conditions. Let 
Qi, bi € {i — 1,2k — i} for i = 1, . . . , fc be fixed, and let the function space Ji be defined by 



This identity also implies that ^ is a positive semidefinite operator on Ji. A condition 
ensuring ^ to be invertible is given in Section 3.1. In the affirmative case the inverse 
operator is given by a so-called Green's function G{t,s) via ^~^f{t) = ^^Q{t,s)f{s)ds. 
Since ^ is positive definite it follows that G{t,s) is positive definite. In particular, the 
matrix defined by 



is positive definite and may be used as the variance of the serially correlated effects Xm ■ 
This specification establishes a link between the covariance matrix a^Ro of the discretizcd 
readings Xm in the model equation (1) and the penalized likelihood equation (2). 

The proposed methodology can be slightly generalized taking .if as the sum of squares 
^f=i of operators measuring different aspects of roughness. The operator may 

be interpreted as a precision and used in the parameterization of a statistical model. This 
is by contrast with standard software for mixed-effects models such as the nlme-package 
[10] in R or the MIXED procedure in SAS, where the parameterization is done in terms 
of variances. In [14] a similar approach was taken for the analysis of longitudinal data, 
and further references may be found in [7], Chapter 8.4. 

The remainder of this paper is organized as follows. Section 2 reviews inference tech- 
niques for the model equation (1). In particular, we present the matrix formulae that 
will be approximated by their operator equivalents. Section 3 provides the mathematical 
contributions of the paper. In this section the operator approximation is introduced and 
refined for the case of equidistant observations, that is, i„ = a-V ^'^^ (& — a) . In particular, 
we derive semi-explicit and numerically stable formulae for the needed computations in 
the case of equidistant observations. In Section 4 the operator approximation is applied 
on the matrix formulae from Section 2. This leads to concrete algorithms that have been 
implemented in an R-package named fdaMixed [8]. 

2. Inference in the mixed-effects model 

This section reviews estimation and inference techniques for the model equation (1). Since 
the derivations of the matrix formulae stated below are standard (see, e.g., [1-3, 13]), no 
proofs will be given. The dimensions are given by 





n.m=\,...,N € IR" 
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where iVtotai = N * M. Based on the covariance matrices G and i? = i?o ® Im we define 
the matrices Ao =1^ + Rq, A = IaTj^j^, + R = Im, and 

Cu = {G-^ + Z^A''^zy^, Cr = A-~^-A-^ZCuZ^A-^, ^ {T^CrT)-^. 

The matrix formulae will be stated such that for moderately sized p and q the computa- 
tional obstacle of their practical implementation lies in the initialization and inversion of 
the A'^-dimensional matrix Aq. The circumvention of this obstacle is the topic of Section 3. 

For known variance parameters cr^, G, Rq, the best unbiased estimate for the fixed 
effects is given by the maximum likelihood estimate 

^ = CpT^CrV = C^r^iA-'y - A-'ZCuZ^A-'y). (4) 

The best linear unbiased predictions (BLUPs) for the random effects u and the serially 
correlated effects x = {xm}m=i,....M arc given by the conditional means 

E[u\y]=CuZ^A~Hy-Tp), E[x\y] = RA-\y -T$ - ZE[u\y]). (5) 

It is generally agreed (sec, e.g., [1] and [7], Chapter 5.3) that the variance parameters 
may be estimated as the maximizers of the restricted likelihood. One of the factors in 
the likelihood is the determinant of Aq =1^ + Rq. To derive the operator approximation 
of this factor we use the representation 

logdetAo= / a„logdet(I + wi?o)dw= / V eJ(i>lAr + i?-^)-iej dw, (6) 
Jo Jo 

where the vectors Cj = {lj=ra}n=i,...,jv G for j = 1, . . . , N constitute an orthonormal 
basis for R^. Using this representation and introducing the conditional residuals r = 
y — r/3 — ZE[w|j/] — E[a::|?;]. the double negative log restricted likelihood is given by 

(2 A^totai - 2p) log <j + M V ej {vIn + Ro^)'^ dv 
Jo ,^1 

+ logdct(Ig + Z^yl"^ZG) -logdetC^ (7) 

+ a-^{r'^r + E[u\y]'^ G-'^E[u\y] + E[x\y]^ R'^'^E[x\y]), 

where it should be kept in mind that C^, r, E[u|?/], E[a;|?/] depend on G and Rq. The 
profile estimate for the error variance has an explicit form, 

= — ^ (r^r + E[u\y]'^G-'nu\y] + E[x\y]'^ R~'E[x\y]) . 

iVtotal - P 

We conclude this section by reviewing some theoretical results on the inference tech- 
niques described above. The errors (3 — [5, E[m|?/] — u, 'Ei[x\y\ — x follow a joint Gaussian 
distribution, and their joint covariance may be derived using [2], Section 2.4. Kackar 
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and Harville [4] show that if the estimators for the variance parameters are translation- 
invariant and even functions of y, then /?, E[u|?/], E[a;|2/] remain unbiased when the esti- 
mates are inserted in place of the unknown variance parameters. As explained by Welham 
and Thompson [16] inference on /3 may be done as x^-tests on twice the log ratio between 
the maximum restricted likelihoods, where the design matrix under the null hypothesis 
is used in the definition of the restricted likelihood under the model. Simulation studies 
done by Morrell [9] suggest that inference on the variance parameters may be done as 
X^-tests on twice the log ratio between the maximum restricted likelihoods, but here the 
formal asymptotic theory appears to be less developed. 



3. Functional embedding of discrete data 

Functional data consist of observations of continuous curves at discrete sample points. 
As an alternative to computations based on spline representations and sparse matrix 
computations we embed the discrete observations into the continuous setting and ap- 
proximate the matrix computations by their operator counterparts. In order to describe 
this operator approximation we first introduce some notation. 

By a discretization of size N of the time interval [a, b] we mean a set of points T = 
{ti, . . . ,tAr} with a < ti < ■ ■ ■ < tN < b. Such a discretization is said to be equidistant 
if t„ = a + ^2jv^ ^ '^)' ^^"^ ™ ^^^^ '■^^^^ associate the mesh length given by A = 
{b — a)/N. To ease notation we implicitly adjoin the points to ~ a and tAr-|_i ^ b to any 
discretization of size A^. 

Given a vector z ~ {zn}n=i,....N G we denote by S'^ the piecewise linear embedding 
of z £ into C([a,6];M), that is, the function that is linear on the segments [i„,f„_|_i] 
for n = 0, . . . , N with (^'z(a) = zi, Sz{b) = zn and Sz{tn) = Zn for n = 1,. . .,N. We also 
introduce the multiplication operator .y£'j- on C([a,&];IR) defined by 

^r./(t) = 4.(i)./W, /eC([a,fe];R), 
where /i = {fJ.n}n=i n G is given from the discretization T via 

{2/(t2+ti-2a), forn^l, 
2/(Wi-t„-i), forn = 2,...,iV-l, (8) 

2/{2b-tN -tN-i), iovn^N. 

In particular, if T is equidistant, then = A^^I. 

Proposition 1. Let a discretization T of the interval [a,b], t € [a,b] and Q G C([a,6] x 
[a,6];M) be given. Assume thatQ{t,-) is twice differentiable on the segments [t„,t„+i] with 
continuous derivatives Q^''\t,). Forz<ER^ there exists {trntn+i) forn^O,...,N 
and Ci G (a,ti), G {tN,b) such that J2n=iS{'titn)zn — J^^ Git, s)S'^{s)S'z{s) ds equals 

2 ^ '(J',Cl)lJ'lZl ^ '{tXN)lJ.NZN 



N . _ 

n— 



where fj. is given by equation (8). 

Proof. The trapezoidal rule of integration [5], Section 7.2, gives intermediate points 
in e {tn,tn+i) such that (?(^, s)^^ (s);?^ (s) ds equals 

ti — a 



2 

n=l 

12 



n=0 

The result follows inserting the piecewise linear functions S"^ and , the first-order Taylor 
expansions at some intermidiate points Ci G (a,ti), ^ {tN,b), 

^-^Q{t, b)S,{h)S, ib) = ^—^Q{t, tN)tlNZN + ^^-^l^g(l) (t, CN)fiNZN, 

by expanding the second-order derivative and by rearranging the terms. □ 

Corollary 1. // the discretization T is equidistant, then there exists ^„ G (^„_i,^„) C 
(tn-i,tn+i) for n = \, . . . ,N such that the approximation error equation (9) equals 

b-a ^ 

~ ~ ~ £,n-l))Q^'^\t,tn)Zn 

— a 



J2 (tn+l - C„) (^^'^ (t, U - G^^^ {t, t^))Zn 
n=l 
N 

67V E(?"-^"-i)(^^''(*'^»)-^'"(^'*"))^" (10) 



12iV 
b — a 



N 



^ " Y,{^n-l-tn-l){Q^^\t,^n-l) - G^^Ht^tn)) Zr, 



UN 

n=l 
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Proof. Equidistant spacing implies that the factors ^„ = N/{b — a) defined in equation 
(8) are constant, and the approximation error equation (9) reduces to 

gjy2 ^ '{t,Cl)l^lZl >{t,CN)l^NZN 

n=0 ^ 

The last sum equals 

n=0 

N ,_ 

+E^(^^'^(*'^"-i)-^''^(*'^"))-- 

n— 1 

By Taylor's theorem there exists ^„ € (Cn-ii^n) such that this equals 

n=l 

The corollary follows centering the terms G'^^^t, •) around ^^^•'(t,t„). □ 
If the matrix D G M^^^ and the integral operator on C{[a, b];M.) are defined by D - 

(^ni ^m)}n,m— 1,..., 

the approximation 

Dz ^{^J^rS A^n)]n=i jveM^, zeM^. (ii) 



{^/(in,im)}n,m=i,...,Af and ^/(t) = J^'' ^(f, s)/(s) ds, then the preceding results suggest 



Green's functions usually possess sufficient smoothness for Proposition 1 to apply 
(see, e.g., [15]), and hence the approximation error in equation (11) vanishes as 
max„=o,...,JV |in+i — in| goes to zero. In case of equidistant discretizations this property 
is refined in Corollary 1. The first term in equation (10) is of size 0(7V~^) and the other 
terms are of size 0{N~^). Perhaps the first term can be used to derive and correct a bias 
arising from the proposed operator approximation, but we will leave this to be studied 
in future work. 



3.1. Explicit operator computations 

To motivate the derivations done in this section we may consider the model equation (1) 
without the fixed and the random effects, that is, y = x + s. In this case equation (11) 
implies the approximation of the prediction equation (5) of the mth serially correlated 
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effect given by 

and the approximation of tlic logaritlimic determinant equation (6) given by 

.1 N 

logdet^o= / V'eJ(t;lAr + i?o ^)"^ej dw 



(12) 



J a 



where Qy is the Green's function for vl + ^^^^ . As shown in Section 4 the matrix 
formulae used for inference in the fuU mixed-effects model equation (1) may be simi- 
larly approximated. In order to develop our computational methodology we derive semi- 
explicit and numerically stable inversion formulae for differential operators of the type 
.5f* = \ + ^ . If the discretization T is equidistant with mesh length A, and the dif- 
ferential operator ^ has constant coefficients, then =1-1- may be inverted using 
Theorem 1 stated below. Boundary conditions play an essential role in this theorem, and 
the reader may want to refresh the definition of the space H. given in equation (3). 



Theorem 1. Consider a differential operator on % given by 

^2fex2/c^ wit/i J, ^C^j^^j ^ be the Jordan canon 



a2k-it^'- ^'{t) + ---+aie'-^Ht)+aoe{t) (13) 



with a2k 7^ 0. Let J ~ diag( Ji, . . . , Jp) G 
ical form of the companion matrix 



( 



C = 






a2k 



1 



a2k 






«2fc-l 

a2k 



l)2fex2fc 



(14) 



Let M e C' 

■>lxfe 



2fcx2fe 



be a non-trivial solution of the matrix equation CM — MJ, and let 
Mij € C^^''^ be the decomposition of the first row of M along the Jordan blocks Jj . Let 



i)eRi^'=, 1^2 = (0 

Fa = {lj_l=aJi 



1)' e 



p2fexl 



and let Fa,Fb G . 



i = l, 
i = l, 



,k ' 
,2k 



Fb = {lj-l=6jj = 1,.. 



T>kx2k 



,k ■ 
,2k 



be given by 



Let Vl = {vi Vl) e 



Fa^ 



Fa 
Okx2k 



&i^2fc^ and let Fa,Fb,We IR2fex2fc defined by 

( 



0fex2/c 

Fb 



W 



Mil 
Mil J I 



MipJp 



MipJl^-^J 
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If the matrix H = FaW exp{aJ) + FhW exp(bJ) is invertihle, then is invertible. In the 
affirmative case the inverse operator is an integral operator fit) = jl^Q^,{t,s)f{s)ds, 
where the Green's function is given by 



G*{t,s) - 



a~^vi e^-p(tJ)H-^ FaW ex-p{{a - s)J)W-'^V2, 



for s <t, 



a-,^viexp{tJ)H-^ FbW exp{{b - s)J)W'^^V2, fort<s. 



(15) 



Proof. The proof follows specializing and condensing [15], Theorem 3. The signs of [15], 
equation (3.15), equation (3.24), should be changed due to a mistake of sign in [15], 
equation (3.9). Wc allow for leading coefficient a2k 1 and have interchanged the indices 
k and p to align with the notation used in the present paper. □ 

Formula (15) is explicit and most satisfactory from a theoretical point of view. But from 
a practical point of view the formula can be numerically unstable since the exponentials 
exp(tJ), exp((a — s)J) and exp((6 — s)J) are weighted against similar exponentials in 
the definition of the matrix H . Imposing symmetry of the Jordan matrix it is, however, 
possible to remove the potential numerical instabilities. 



Proposition 2. Suppose that the characteristic polynomial 

a2kz'^^ + a2k-iz^^^^ H h aiz + qo = 



(16) 



for the differential operator (13) has 2k distinct roots rj^ ,rj^ , . . . ,7]^^ ,ri^ G C such that 
the real values of the k eigenvalues r^^ , . . . ,rj^ are non-positive and the real values of 
the k eigenvalues r]^,...,T]'^ are non-negative. Then the Jordan canonical form of the 
companion matrix equation (14) is diagonal with block diagonals consisting of eigenvalues 
with non-positive and non-negative real values, respectively, 



J- 
Ofexfe 



Ofexfe 

J- 



and the matrix W : 



\w+) 



J_ =diag(77i ,...,11^ ), J+ =diag(?7+,...,7y+), 
■^2kx2k decomposed via W-,W+ G £2kxk ^g^j^g^ 



1 



-\2k~l 



1 

Vk 

iVk?"-' 



Furthermore, define vi = (1 ■ ■ ■ 1) e 



p 1 X fe 



W+ 



V2 



1 

rit 



+ \2k~l 



(0 ••• 1)^ e 



1 



p2fcxl 



via W ^V2 



^2k 



I, and the vectors (l)^{t),ip^{t) £ 



p 1 X fe 



for t G [a, b] and G No by 



pfcx 1 



{viJf-viJ^lc-^^-'^^+iFbW+y^FbW-C^^-'^-^-) 



■ {Ikxk-e^'^''^'-{FaW.r^FaW+e-^^~''^'+{FbW+r^FbW.e^''-''>'~)-^ 
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and 

Then the fith partial derivative dfG^{t,s) of the Green's function defined in equation (15) 
may be rewritten as the numerically stable expression 



(17) 



/.^(Oe(*-"'^-(v_+e("-"'^-(FaM^_)-iFaM^+e-(^-")'^+w+), fors<t, 
~iP,,{t)c-^-'-'^-^+{v++c-^''-'^-'+iFbW+)-^FbW-C^''-'^-'-v^), fort<s. 

Proof. From equation (15) we have that dfG*{t,s) equals 

a-^viJ^" eMtJ)H^^ FaW eyipiia - s)J)W~\2, for s<t, 
-a-^viJ''exp{tJ)H-^FbWeyi-p{{b~ s)J)W-^V2, for t < s. 

The crux of the reformulation of this representation lies in the inversion of the matrix 
H = FaW exp{aJ) + FbW exp{bJ). To this end, we write He~*'^ and e^'^H^^ as block 
matrices with k x fc-blocks, 



-[FbW.e('-'^J- FbW+e('>-'^J+ ) ' ° ^ " l^A^i A2 



Using elementary matrix algebra we find that wiAn + W2A21 for general wi, W2 G M^^*^ 
equals 



iwi-W2e-'^''-*^-^^iFbW+)-^FbW-e^''-''^^-) 
{FaW-e^"'*''^- - FaW+e-'^''-'''>^+ {FbW+)-^FbW-e^''-'^-^'y^ . 
Inserting this above we have that d^Q^:{t, s) for s <t equals 

-If T^\f^il Ai2\ /FaW/^cf"^-^)^- FaW+C^''-''>-^+\ f 



(18) 



which equals 

a^^{virAn+viJ'lA2,){FaW-e^''-''>'-v.+FaW+c^''~'^''^v+). (19) 

Combining equations (18) and (19) and rearranging the exponential factors we arrive at 
equation (17) for s <t. The reformulation is done similarly for < < s. □ 

Remark. From the viewpoint of statistical modeling, the results in [15] are more general 
in two valuable ways. Firstly, the boundary conditions separately given at the end-points 
of the sample interval via the matrices Fa, Fb in Theorem 1 may be given in form of 
linear combinations of the curve and its derivatives at a and b via general Fa and Fb- In 
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particular, boundary conditions enforcing periodicity may be stated. But to derive the 
numerically stable formulae stated in Proposition 2, we have refrained from this possibil- 
ity. Secondly, the results in [15] are given for matrix- valued functions. This generalization 
allows our methods to be extended to multivariate functional data analysis. 

In the following theorem the explicit inversion formula is applied to derive a simulta- 
neous computation of i9f (I -f- A^)~^ S'z{tn) for n — 1, . . . , N , where z G M^, that easily 
may be implemented with computational complexity 0{N). Furthermore, the inner inte- 
gral in the approximation equation (12) of the logarithmic determinant may be explicitly 
computed for Lcbesgue almost all v G [0, 1]. In the statement of the theorem we denote by 
the element-wise multiplication of matrices or vectors of the same dimension. Unless 
specified otherwise the operation is performed after ordinary matrix multiplications. 

Theorem 2. Suppose the discretization T is equidistant with mesh length A = (b ^ 
a)/N, and assume that the operator in equation (13) given by =I + A^ satisfies the 
conditions of Proposition 2. Denote by the Green's function for let J_, J-|_, W^, 
W+ , v^, , (j)^ (t) , (t) be as defined in Proposition 2, and let ■J- 1 , ^- j C-i- 1 C+ ; C+ S 
R*^^^ be defined by 

rexp(A7£/22-11 ^ f l-cxp(-A7?+/2) | 

,0 ^ f l-(l-A7y,r)cxp(A7?~) ] ,0 _^ f exp(-A7?+)-l + A7?+ ) 

I A(,,r)2 4+ \ Ai4r i.=i,...,.' 

1 ^ f exp(A77,r)-l-Ar;-) ] ,i ^ f 1 - (1 + Ar;+) exp(-A77+) ] 

For z = {zj}j=i,....Ar G M.^ the fith derivative 9f (I + A^)^^(^z{tn) taken at the sample 
point tn is given by 

n-l 
n 

+ 0,.(i„) ^e(*"-*^)^- («_ (l.^iC- + lj>id))^, 
i=i 

TV 

- VM(i„) 5]e-(*-*")^+K {l,<Ne+ + 1,=JV?+))^, 

j=n 

N 

-V'M(in)l„<iv e-(*^-^-*")'^+(w+0Ci)^, 

j=n+l 

n-l 

+ 0p(t„)e(*"-")-'-(^^aW^_)-ii^aT^+l„>i ^e-(*^-'^)-^+(z;+ 0e^)z, 
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+ 0M(in)c^*"""^^-(^aM^-)-'i^aM^+^C-(*-i-")'^+(i;+0(lj = lC++lj>iei))Zj 

N 

- ^^(i„)e-(^-*")-^+(i^f,T^+)-ii^f,W^_ ^e(''-*^>i)''-(z;_ (l,<ive° + WC-))^, 



AT 

j=n+l 

Concerning the log determinant assume that the operator in equation (13) given by 
^* = vl + for fixed v £ [0, 1] satisfies the conditions of Proposition 2. Let the matrices 
G]R'=^'= be defined by 



j ' > i,j=l,...,k 



ti 'Ij ) ) i,j=l,...,k 

and let the matrix B e M'^^'^ be defined by 

{FaW-)-^FaW+c-''''-''^''+iFbW+)-^FbW^ 
{Ikxk~c^''-''^-'-{FaW-r'FaW+c-^''~^^'+{FhW+)-'FbW^)-\ 

Denoting by t > the leading coefficient of J£ , then the integral j'' Q^{t,t) dt equals the 
sum of the following 8 terms: 

I ^ Nt'^viv^, 

II = T-'vi{{FaW-y'FaW+ 

/// = -T-\i{{FbW+)-^FbW- A+_)v-, 

IV = -T-iz;i((FbW^+)-iFf,M^_e(^-'^)^- {F^W-Y^FaW^ 

VI = r-^i(Bc(''-"''^- {FaW^Y^FaW+ 
F// ^ -T-^vx(iFbW+Y^FbW-(i^^-'''>-'-BQ Aj^Y'"^. 

VIII = -T~^Vx{{Fy,W+Y^FbVV-4-^-''^-^-BB^^-'''^-^- {FaW-Y^FaW+(z)A++Y+. 
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Proof. Since the characteristic polynomial has distinct roots 771, . . . , 772^, the Jordan 
canonical form of the companion matrix is diagonal, and equation (15) implies that 
d^{I + A^)-^,^,{t) equals 

a-^viJ'' CJip{tJ)H~^ FaW c^p{aJ)I^J e-"'''^^(s) ds • (W^-^W2)i| 



l,...,2fc 



t J i=l....,2fe 



Since the function S'z is piecewise linear, the above integrals can be explicitly evaluated 
over the intervals [tj,tj^i]. For j = 0,N, we have 

/\--^.(.)d. = e--i:i^^P(z^,, 
Ja m 

re--^.(.)d. = e— i^^^Piz^,,, 

and for j = 1 , . . . , — 1 , we have 

(■tj+i 

e-^'''^^(s)ds 

.A 

Jo 

= e-*^"' / e-'''''(l- A-\s)dszj+c-*^''' / c-"-^A-\sdszj+i 
Jo Jo 

-t,„. cxp(-A77,)-l + A7^, -t,„. 1 - (1 + A?/,) cxp(-A77i) 

Arranging the eigenvalues as r]^ , . . . ,7]^ ,t]^ , . . . ,ri^ and inserting the definition of ^_ , 
C+, C+, 61, wc have that (I + Aif )-i(r^(i„) equals 



n-l 



ln>i E a2-.^«i J-e*"^i7-F.W^e(-*^)^ ( ^^^^^^^^ 



N 



1„<. Y: a2-.^^i^^e*-^F-F,T^e(''-*-^)'^(^-®^;^T^ 



i=n+l 
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The exponential factors on the terms may be assimilated in the expo- 

nential factors before the large parenthesis using tj+i — tj = A for j = 1, . . . , iV and 
t2 — ti = t]\r^i — tjv = A/2. Thereafter the terms in these sums are of the same type as 
in equation (19) with w_, replaced by w_ , etc., and the formula for 
d^{I + A^)~^S'z{tn) follows by invoking the same reformulations as used in the proof of 
Proposition 2. 

Finally, we consider the Green's function CJ* for = vl + A^ . The differential oper- 
ator has leading coefficient cx2k = At, and inserting s = ^ in the first part of equation 
(17), we find that G*{t,t) equals 

A-^T-\vie'^''-''>''- ~vie~'^''-'^-^+{FbW+y^FbW-) 
{h>ck-e'-''-''^'-iFaW.r'FaW+e-^''-^^'^iFtW+)-'FhW.y' 

To remove the possibly exploding exponential factor e"^^"*)''" in the first factor, we 
invoke the matrix formula (/ — X)~^ = I + X{I — X)~-^ on the second factor and rear- 
ranging the exponential factors. Doing this G*{t,t) is rewritten as the numerically stable 
expression 

A-'T-^{vi-vie-^''-'^-^+{FbW+)-^FbW^c^''-*^''-) 
(w_ 4-c(*-")^- {FaW-y^FaW+c~^*~'''^-^+v+) 
+ A-^T-\viC^*-''^-^- -vie-^''-''>'^+{FbW+)-^FbW^e^''~'''>'^-) 

iFaW.)-^FaW+e-'^''~''^''+iFbW+)-^FbW^ 
iIkxk-e'^''-''^'^-iFaW^)-'FaW+e-^''-''^'^+{FbW+)-^FbW.y^ 

This expression is expanded into the sum of 8 terms, which all may be explicitly integrated 
over the interval [a,b]. For instance is the integral over the second term given by 

fa 

A- V-^vie^*-'^)"'- iFaW-)-^FaW+e-^'-''^-^+v+ dt, 

which equals T-^vi{{FaW-)-^FaW+ A_+)v+. □ 

Remark. The predictors E[a;m|j/] rnay be seen as the predictors E[a;^*|2/] for the func- 
tional parameters x^' evaluated at the sample points i„. The formulae stated in Theo- 
rem 2 may be extended to functional representations for E[a;^'|y]. Doing this the predic- 
tions between sample points will be given as linear combinations of exponential functions. 

Remark. If the kernel G{t,s) is constant, say G{t,s) = A, then the operator approxima- 
tion 

^(«I + ^if)- (t,) dv = £ £ -^^1^^"°^ dtdv = log(l + NX) 
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gives the exact log determmant of {l„=m + G(tn,tm)}n,m = Iat + {^}n,m.- The particular 
construction of the embedding operator (o^ was chosen to achieve this property. 

A fundamental difference between our operator methods and the smoothing spline 
technology lies in our dependence on boundary conditions. Whether boundary conditions 
are desirable in statistical modeling depends on the data situation at hand. If wc have 
additional knowledge implying particular boundary conditions, then this may be used 
in the statistical model. However, in many data situations such additional knowledge is 
not available, and the requirement to specify boundary conditions may be disturbing. 
Here our advice is to use Neumann-type conditions. Although the covariance function 
G{t,s) is not defined for Neumann conditions as noted in the following example, this is 
possible due to the regularization induced by the measurement noise; that is, 1 + .^^'^^ 
is non-singular by construction. 

Example. For ~ Xdt wc have Ji' = Jif'^Jlf = —X^df. Consider the following two sets 
of boundary conditions: 



Thus, the Laplace operator with boundary conditions (Bl) leads to the Brownian motion, 
and the Laplace operator with boundary conditions (B2) leads to the Brownian bridge. 
The Laplace operator with Neumann boundary conditions 6'^^-'(a) = = is not 

positive definite. Even so, this operator can be used in a statistical model, where it implies 
an improper prior for the serially correlated effects in terms of a Brownian motion with 
a free level. 

To compute the approximative log likelihood we find the Green's function Qy for vl + 
^jf^^. In case of the Brownian, motion equation (17) gives 



{Bl): 61(a) =6*^1) (6) = 
Wc have if-i6'(i) = ^lg{t,s)6{s)As with 



{B2): 0{a)=e{b) = O. 




for boundary conditions (Bl), 
for boundary conditions (B2). 






In case of the Brownian bridge, equation (17) gives 
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sinh((((t A s) - a)/{Xy^b~~^))y/m) sinh(((& - (t V s))/{Xy/b~^j)y/Nv) 



L 



2A 



1 



In both cases the double integrals jl^Gvit,t) dt dv can be computed giving explicit for- 
mulae for the operator approximation of the matrix determinants. In case of an equidis- 
tantly sampled Brownian motion, we have 



logdet{l„=™ + tj(i„,im)}„,m=i,...,jv ~ log(cosh(A ^y/b-aVN)) 
and in case of an equidistantly sampled Brownian bridge, we have 

l0gdet{ln=m + y[tn,tm,)}n,m=l,...,N ^^Og' 



4. Approximative inference 

In this section we combine the matrix formulae listed in Section 2 with the operator 
approximation developed in Section 3. The obstacle in the matrix computations is the 
inversion of the matrix Aq = Iat + i?o G M^^^. Here i?o = {Q{tn,tm)}n,m=i....,N is defined 
via a discretization T = {ti , . . . , ijv} and the Green's function Q for a differential operator 

-1^1=1^1 ^l- 

The maximum likelihood estimator and the BLUPs given in equations (4) and (5) are 
approximated using the block structure A = Aq® Im, the identity A'f^^ z — z — R^A'f^^ z 
for z S and the approximation 

RqA^^z = A^^Roz = (I + RoY'z « {(I + '^)- V.(t„)}„=i,...,iv. 

Note that this approximation is applied both on the individual sample vectors e 
and on the sections of the columns of the design matrices F and Z . The approximation 
of the logarithmic determinant equation (6) in the restricted likelihood equation (7) has 
already been stated in equation (12), and the quadratic form of the serially correlated 
effects is approximated by 

M N 

E[x|j/]^i?-iE[x|y] « ^ ^E[a;™(t„)|y]^.^^iifE[x™(t„)|y]. 

m— 1 n— 1 

Furthermore, for an equidistant discretization with mesh length A, we have 

L M N 

E[x\y]'^ R-'E[x\y] ~ ^Y.Y. E(=^EK„(t„)|y])T(J^iE[x„(t„)|2/]). 

m— 1 n—1 
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If the discretization T is equidistant, then semi-expUcit and numerically stable formulae 
for the above approximations are given in Section 3.1. For general discretizations the 
operator approximations may be found as numerical solutions to ordinary differential 
equations; for example, the function / = (I + ^■^^^)~^ S'z € H obeys to the differential 
equation / + ^^^J^f = S^. 
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